% function demo_cavity
clear all; 
global V T E TE ET posV posT posE;
caseNum = 4;  % for driven cavity
tic;
% for chebyshev grid:
%%%%%%%%%%%%%%%%%%%%%%%%%  mesh 1  %%%%%%%%%%%%%%%%%%%%%%%%%%
% [V,T] = chebytri(12);  
%%%%%%%%%%%%%%%%%%%%%%%%%  mesh 2   %%%%%%%%%%%%%%%%%%%%%%%%%%
% V = [0 0; 1 0 ; 1 1; 0 1; .5 .5];
% T = [1 2 5; 2 3 5; 3 4 5; 4 1 5];
% [V,T] = refine_uniform(V,T); 
% [V,T] = refine_uniform(V,T);
% [V,T] = refine_uniform(V,T); 
%%%%%%%%%%%%%%%%%%%%%%%%%  mesh 3   %%%%%%%%%%%%%%%%%%%%%%%%%%
% % for another test mesh in unit squre;
% [p,e,t] = initmesh('squareg','Hmax',0.45,'Hgrad',1.99);
% posV = size(p,2); posT = size(t,2);
% V(1:posV,:) = (p(1:2,:)'+1)/2; 
% T(1:posT,:) = t(1:3,:)';
%%%%%%%%%%%%%%%%%%%%%%%%  mesh 4 %%%%%%%%%%%%%%%%%%%%%%%%%%
V = [  0 1;    0.2 1;    0.5 1;     0.8 1;     1 1; ...
    0 0.98; 0.2 0.98; 0.5 0.98;  0.8 0.98;  1 0.98; ...
     0 0.5;  0.2 0.5;  0.5 0.5;   0.8 0.5;   1 0.5; ...
     0 0.2;  0.2 0.2;  0.5 0.2;   0.8 0.2;   1 0.2; ...
       0 0;    0.2 0;    0.5 0;     0.8 0;     1 0  ];
T = [ 1 6 2; 2 6 7; 2 7 3; 3 7 8; 3 8 9; 3 9 4 ;4 9 10;4 10 5; ...
    6 11 7; 7 11 12; 7 12 8; 8 12 13; 8 13 14; 8 14 9; 9 14 15; 9 15 10; ...
    11 16 17; 11 17 12; 12 17 18; 12 18 13; 13 18 14; 14 18 19; 14 19 15; 15 19 20; ...
    16 21 22; 16 22 17; 17 22 23; 17 23 18; 18 23 19; 19 23 24; 19 24 20; 20 24 25];
[V,T] = refine_uniform(V,T);%[V,T] = refine_uniform(V,T);

posV = size(V,1); posT = size(T,1);
areas = tdata;
d = 5*ones(posT,1); % d(10)=10; % the last sencense if for test;
[cr,desc,asce] = get_auxillary_mat(15);
[Mat1,Mat2,Mat3,i_pattern,j_pattern] = ns_common_mat(5,6,desc);

% [dofs,n_dof] = sort_dof(d,T,E,TE,ET,posV,posT,posE,0,cr);
% [H,row_idx] = smooth_C1(dofs,V,T,E,posE,ET,TE,d,cr);
[dofs,n_dof] = sort_dof_dis2(T,posT,d);
% [H,row_idx] = smooth_C1_dis(dofs,V,T,E,posE,ET,TE,d,cr);
H = cr_matrix(dofs,V,T,E,posE,TE,ET,desc,asce,cr,d,1);

% linear part of system matrix
K = bending_bb_2(dofs,V,T,posT,d,areas,desc); %check it
M = mass_bb_2(dofs,V,T,posT,d,areas);
% boundary condition .
bdr = sort_border(V,E,ET,posV,posE);
[B,G]= inflow_bc(V,T,TE,ET,dofs,d,cr,bdr,'g1','g2',caseNum);
% [B,G]= get_navstk_bc(V,T,E,TE,ET,posV,posT,posE,dofs,d,cr,'g1','g2',caseNum);
b = bnet_2(dofs,V,T,posT,d,'h',caseNum);
% set the parameters
Max_linear_iter = 3; linearTol = 1e-8; epsilon = 1e-6;
Max_newton_iter = 12; newtonTol = 1e-8; 
renold0 = 100;renold_end = 5000;

p = size(H,1); ZERO = zeros(p,1);
L = [H;B]; G = [ZERO;G];
% then solve the ns problem using init renold number.
renold = renold0;
% get the solution of stokes problem for initial value;
c0 = lagrange22(K/renold,M*b,L,G,epsilon,Max_linear_iter,linearTol);
% and solve ns using c0;
[c,cvg,diff] = navstk_static(V,T,posT,dofs,d,K,M,L,G,b,Mat1,Mat2,Mat3,i_pattern,j_pattern,c0,renold,epsilon,linearTol,Max_linear_iter,newtonTol,Max_newton_iter);
iter_end = 0;
c_good = c;  % save for plot;
renolds = zeros(100,1); % the stack simulator
renolds(1) = renold_end; pos = 1;
while ~iter_end  % iter until some condition meet.
    if cvg
        fprintf('----For Renold = %f-------Convergenced-------------\n',renold); 
        % save the better result for plot;
        c_good = c; 
        Renold_good = renold;
        if renold < renold_end % if have not get the desired renold,              
            c0 = c;       % init value     
            % save the current convergenced renold as renold0
            renold0 = renold;  
            % get the next renold from vector renolds
            renold = renolds(pos);
            % the position of index reduced itself
            pos = pos - 1;
            [c,cvg] = navstk_static(V,T,posT,dofs,d,K,M,L,G,b,Mat1,Mat2,Mat3,i_pattern,j_pattern,c0,renold,epsilon,linearTol,Max_linear_iter,newtonTol,Max_newton_iter);
        else   % if the disired renold if meet ,plot the solution and quit now.
            iter_end = 1;
        end
    else
        fprintf('----For Renold = %f-------Not Convergenced---------\n',renold);
        if renold-renold0 <= 10  % actually, cant's comput further more
            iter_end = 1;
            fprintf('****************************************************************\n');
            fprintf('Sorry!I can not compute the cavity solution of renold number bigger than %f.\n',Renold_good);
            fprintf('****************************************************************\n');
        else
            c0 = c_good;
%             first store the current renold for later use
            pos = pos + 1;
            renolds(pos) = renold;
%             and half it
            renold = (renold0 + renold)/2; % pick the mid-point of current unconvergence renold number and the know convergenced renold number
            [c,cvg] = navstk_static(V,T,posT,dofs,d,K,M,L,G,b,Mat1,Mat2,Mat3,i_pattern,j_pattern,c0,renold,epsilon,linearTol,Max_linear_iter,newtonTol,Max_newton_iter);
        end
    end
end
% plot the computational mesh.
% figure;plot_t(V,T(1:posT,:));title('The computational mesh.');
plot_cavity_stream;
plot_cavity_votocity;
toc;